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(54) TiUe: IMPROVED METHODS FOR PERFORMING RESERVOIR SIMULATION 
(57) Abstract 



A method for performing reservoir simulation by 
solving a mixed implicil-IMPES matrix (MIIM) equation. 
The variable implicit reservoir model comprises a plurality 
of cells including both implicit cells and IMPES cells. 
The MIIM equation includes a firet scalar IMPES equation 
for each of the IMPES cells and a set of implicit 
equations for each of the implicit cells. The reservoir 
simulation method comprises the following operations: (a) 
constructing a global IMPES pressure equation fiXMn the 
MIIM equation; (b) solving the global IMPES pressure 
equation for pressure changes; (c) computing first residuals 
at Che implicit cells in response to the pressure changes; 
(d) detemiining improved saturations by solving the total 
velocity sequential equations at the implicit cells using 
the first residuals; (e) computing second residuals at the 
implicit cells and at a subset of the IMPES cells that are 
in flow communication with any of the implicit cells in 
response to the improved saturations; (f) determining if 
a convergence condition based on the second residuals is 
satisfied; and (g) repeating steps (b) through (0 until the 
convergence condition is satisfied. After the convergence 
condition is satisfied, a final solution estimate for the 
MIIM equation is computed from the pressures changes 
and the improved saturations. The final solution estimate 
may be applied to detemiine the behavior of the reservoir 
model at a future discrete time value. Alternative to stq) 
(d), improved saturations and improved pressures may be 
computed by peifonning one or more iterations with a 
selected preconditioner at the implicit cells. 
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TITLE: IMPROVED iMETHODS FOR PERFORMING RESERVOIR SIMULATION 



Field of the Invention 

The present invention relates to reservoir siomlation, and in particnlar, to methodologies for performing 
reservoir simuktion by solving an iirq)licit matrix equation or an uxq)licit-IMPES matrix equation. 

Background of the Invention 

In an attempt to understand and predict die physical behavior of reservoirs (such as petroleum reser* 
voirs), reservoir engineers and scientists have generated various mathematical descriptions of reservoirs and the 
fluids they contaixL These mathematical descriptions are often expressed as coupled sets of differential equa- 
tions. Smce it is quite often inq)ossible to obtain solutions of the differential equations in all but the simple 
cases, the differential equations are discretized in space and time, and the resulting difference equations are 
solved using various immerical simulation techniques. For exaiiq)Ie, the following difference equations repre- 
sent the volumetric accumulation of oil and water in a particular cell (i.e. cell t) over the course of a timestep 
from time index n to n+1 assuming rock and fluid incompressibility in a one-dimensional reservoir 



(Bl) 



(B2) 



where ^ is die timestep size; 
y 

' is the volume of cell i; 
^ is porosity, i.e. pore volume per cell volume; 

^^o)i |3 ^ saturation of oil at cell i, Le. the fraction of the pore volume occupied by oil in cell i; 
(Sw)i ^ ^ sanuation of water at cell i, Le. the fraction of the pore volume occiqpied by water in cell 

i; 

B B 

^ and * are the formation vohmie factors (FVF) for oil and water respectively; 

(/^o)/-! ^ (Po)i ^ iPo)M are oil pressures at cell i-1, cell i, and cell i+1 respectively; 

iPw)tA ^ iPw)i ^ (Pw)m water pressures at cell i-1, cell i, and cell i+1 respectively; 

(^o )i is the rate of oil injection into cell i, and takes die value zero at most cells and takes a negative 
value at cells which interact widi a depletion well; 

(9«r)i [3 Hie late of water injection into cell i, and typically takes a zero value except at cells which 
interact with an injection or depletion well; 
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and W*** represent a quantity x evaluated at time indices n and n+1 respectively, where the 
former is known information, having been determined from previous computations, and the later is an unknown 
to.be solved for by some con^tational method; and 

and ^^^^ represent quantities which are to be cvahiated at time index n or n+1 subject to user 

selectiorL 

The oil transmissibility-mobiUty factors ^"^'^'^^i and ^'^"^'->^ are defined as 



(B3) 



(B4) 



where A is the area normal to the axis of the one-dimensional reservoir; 

(^fl)/+>$ is the mobility of oil in transit between cell i and cell i+1 ; 
(^0 )/->^ mobility of oil in transit between cell i and cell i- 1 ; 



is the position of the kth cell along the one-dimensional axis. 

Similar dcfinitioiis apply for the water transmissibility-mobility products and Tne 

difference equations (Bl) and (B2) above arc augmented widi several auxiliary relations as follows: 

Pw-Po =^c('S'o)^ (B6) 

^w-^w(Pw*^w)^ (B8) 

p 

Relation (B5) follows from the definition of saturation. Capillary pressure * which is defmed as &e difference 

M 

20 in pressure between water and oil is a known function of oU saturation. Oil mobility * is a known function 

of oil pressure and oil saturation. Water mobility is a known function of water pressure and water sattira- 
tioiL 

Smce oil mobility is a function of oil pressure and oil saturation " , and Aese later variables 
are defined at ceil centers, a question arises as to the proper means of evaluatiiig the in-transit oil mobilities 

25 (^«^f+)i 3jjj| /According to the midpoint weighting scheme, the in-transit oil mobility is defined to 

be the average of the mobilities at the two afiected cells. For example, 
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15 



20 



where is evaluated using the oil saturation and oO pressure ^^'^^^ prevailing at cell i, and 

(Mo)m is evaluated using the oil saturation ^^•^'+* and oil pressure ^^«^'+* prevailing at cell i+L Alterna- 
tively, according to die upstream weighting scheme, the in-transit oil mobility may be defined as the oil mobility 
at die upstream cell of die two affected cells, where the upstream cell is defined as die cell with higher pressure 
(since fluids flow from high pressure to low pressure). For exan^le. 



^ '^'*^~\{MXi> otherwise. 



(BIO) 

If du pressure variables and transmissibility-mobili^ factors in Equations (Bl) and (B2) are evaluated 
at Ac new time index, ie. ^~^~''"'"^,Equations(Bl)and(B2)takefliefiarm 



(^.)r=;gjkr-(al 



(BID 



(B12) 

The transmissibility-mobiHty factors and the phase injection rates are functions of saturation and pressure, and 
are evaluated at die new time level n+1. Hius, Equations (Bl 1) and (B12} are non-linear in die unknown vari- 
ables 

iPo)f (j>o):^ 

(pj)ti cPw)r' iPwYiti 

and^^"-^'-'', (B13) 

Equations (Bl 1) and (BI2) may be expressed in terms of a reduced set of unknown variables using relations 

may 



(B5) and (B6). For example, die variable ^^"^^ may be replaced by ^ ^^^^^ . Similarly, ^^"^^ 



be replaced by ""^^ K Unis, Equations (Bl I) and (B12) may be expressed in temas of die 

following reduced set of unknown variables: 

Assuming that there are N cells in the reservoir being modeled, Equations (Bl 1) and (Bi2) describe a coupled 
non-linear system of 2N equations (two equations per cell) widi 2N unknowns - each cell contributes an un- 



3 
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known pressure ^' and an unknown saturation ^ * .An iterative method such as Newton's method 
is generally required to solve such systems. 

Let vector <^ be tiiie vector of 2N unknowns for the system. Define a set of 2N functions ^, 

j = 1, 2, 3, 2^ ^ functions per cell, as follows. A first function ^"^^^^ for cell i is defined by the 
5 expression which follows from subtracting the right-hand side of Equation (Bll) from the left-hand side of 

Equation (Bl 1). A second function ^ for cell i is defined by the expression which follows from sub- 

tracting the right-hand side of Equation (B12) from the left-hand side of Equation (B12). Let J * 

be the corresponding vector function whose component functions are the functions . The system given by 
Equations (B 1 1) and (B12> may be equivalently expressed by the equation 

10 = (BIS) 

It, the solution X^X* of the system given by Equations (Bll) and (B12) conresponds to the zero of Equa- 
tion (B15). Equation (B15) may be referred to as a fiiUy in^licit equation or a nonlinear implicit equation since 
none of the unknowns (B 14) may be explicitly computed from known data. Hius, any method of solving equa- 
tion (B15) may be referred to as a fully implicit me^d. 
15 Newton*s method prescribes an iterative method for obtaining the solution of Equation (B15). Given a 

X f 
current estimate * of the solution, the function is linearized in the vicinity of this current estimate: 

Y^df{X,){X-X,)^f{X,)^ 

where ^^^*^ represents the Jacobian matrix of flie fimction ^ evaluated at ~ , and ^^^^^ repre- 

f X 
sents die evaluation of function ^ at the current estimate. The next estimate of the solution is obtained 

X 

20 by setting vector Y equal to zero and solving for argument X. Thus, the next estimate satisfies the matrix 
equation 

X X 

By solving Equation (B17) for successively increasing vahies of the index k, a sequence of estimates ^ , S 

...^ ... is obtained which converge to the solution of the noiilinear system (BIS). 
25 Equatbn (BIT) is referred to herein as an implicit matrix equation. A linear equation solver is used to 

solve the unplicit matrix equation (B17). The right-hand side vector ^ " ^* ^ and the Jacobian 

matrix ^^^^^ are supplied as ix^ut data to the linear solver. The linear solver returns the solution vector 
y 

M ^ inq)licit matrix equation (BIT). The computational effort of a Newton*s mediod sohition of die 
nonlinear implicit equation (B15) depends on (a) the number of Newton iterations to achieve convergence of die 
X 

30 sequence ^ , (b) the average corrqnitational effort expended by the linear solver to solve the in:9>licit matrix 
equation (BIT), and (c) the computatioimi effort required to i^date tihe matrix equation as improved solutions 

4 
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are obtained. While most of the computational effort per Newton iteration is associated widi solving the matrix 
equation* tfie effort required to update die matrix equation is also significant Thus» any in^rovement in the 
computational efficiency of die linear equation solver will have a coiiespondtng effect on die efficiency of the 
Newton's mediod sohition. 

As described above, the nonlinear implicit equation (B15) arises from the choices ^"fi^'^'^^ jn 

Equations (Bl) and (B2) above. Another plausible set of choices is given by CC^n-^l ^"''.where- 
upon Equations (Bl) and (B2) take die form 



-^(or = -fekr-(al 



. (B18) 

+(^j7=^[(or-(a] 

(B19) 

The saturations and pressures at time-index n comprise known data (having been determined from previous 
computations). Thus, the transmissibility-mobility functions evaluated at time-index n comprise known con- 
stants. Equations (B18) and (B19) are therefore linear in the unknown variables 

iPo)H iPo)r 

t 9 f 

(Pw)m ^ (Pw)t ^ iPw)M ^ (B20) 

One method for sohring die linear system of Equations (B18) and (B19), Le. die so called Implicit- 
Pressure Explicit-Saturation (IMPES) method, is motivated by the following reduction of Equations (B18) and 
(B19). Since the saturation variables obey relation (B5\ the Equations (B18) and (B19) may be combined so as 
to eliminate the unknown saturation variables. In particular* Equation (B18) may be multiplied by the oil for- 

B B 

mation volume factor ° , and Equation (B19) may be multiplied by die water formation volume factor ^. 

The resulting equations may be added together to generate the following linear equation involving only the pres- 
sure unknowns: 

Boi^)i}i[(Po)::: -(p.)n - BMJi^[(Po)r -(pxh 

+i^.(0?+^w(9»)7 =0 . (B21) 

Equation (B21) is refened to herein as an IMPES pressure equation. The capillary pressure relation (B6) may 
be used to eliminate die water pressure unknowns under die assunqition diat capillary pressure does not change 
during the timestep: 

(pjr=(p.)r+w.);], 

where j represents an arbitrary ceil index. When Equation (B21) is written for all N ceils in the reservoir, the 
ensuing system, herein referred to as die IMPES pressure system, has N equations and N unknowns - one un- 
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knownpitssutc^"^'*' per cell Because the IMPES pressure systm b linear 
Hiikii«>wns it nay 1» sohred nuKA fiater to &e ftlfy 

Again 8 linear equation solver nay be invdced to solve the Hie solution 

vector of the IMPES pressure system specifies the pressure at every ceU in the reservoir. Tbe 

5 unknown saturations and in Equations (B18) and (B19) may be dctemiined by substituting 

the pressure soh«ion vahws the left-hand sides of Equations (B18) and (B19). Since the satura- 

(S )**' 

tions and are known ftom previous computations, the unknown saturations • ' and 

(^w)r' a^y be computed expUcitly. Thus, the IMPES procedure involves two steps: a first step in which 
pressures are computed impUcitly as the solution of a linear system; and a second step in which satarations are 

10 conqjttted explicitly based on the pressure sohition. 

The example of a one-dimensional model discussed above represents a greatiy simplified description of 
a complicated physical situation. More realistic modeb invoWe (a) a two-dimensional or three-dimensional 
array of cells, (b) more than wo conserved species, (c) more than two phases, (d) compressible fluids and/or 
rock substrate, (e) non-uniform ceU geometry and spacing, etc. In addition, the difference equations of the res- 

15 ervoir model may not necessarily arise ftom a fluid votome balance. In otiier approaches, difference equations 
niay be obtained by performing, e.g., mass or energy babnces. WWfc pressure is quite often one of the vari- 
ables being sowed for at each cell, the remaining variables need not necessarily be samrations. For example, in 
other formulations, the remaining variables may be mole factions, masses, or other quantities. 

CHvcn a reservoir with M conserved species, a conservation law may be invoked to write a set of M 

20 difference equations describing the physical behavior of each of the conserved species at a generic eeU i. (The 
use of a single index i to denote a generic cell does not necessarily imply that die reservoir model is one- 

dimensionaL) Hie set of equations may generaUy be expressed in tenns of the pressure ^» of some base spe- 
cies (often ott), and (M-1) complementary variables such as saturations, mole fractions, masses, etc. These 
complementary variables will be referred to herein as gencrali2ed sattirations. 

25 The discussion of die fiiUy impUcit metiiod and the IMPES method presented above generalizes to 

more realistic models. TT«M difference equations fi« the genetic ceUigemmdlyinctodefimctions such as mo^ 
biUly. formation volume fector, pore votame. injection rate etc, which depend on pressure and/or the general- 
ized sanirations (Le. complementary variables). The fiilly implicit equations result fiom evaluating such func- 
tions at die new time mdex n+l. The fully impKcit equations are generally non-linear, and dms. requite an it- 

30 erative method such as Newton's method for their solution. 

The IMPES formulation starts from evahiating functions of pressure and/or die conqilementaty vari- 
ables at the oM time index n. ihus. the M difference equations particubrize to a set of linear equations in the 
unknown pressures and unblown genaaKzed saturations. An auxiliary reUtion analogous to relation (B5) may 
be used to combine the set of linear equatioiB into a single equation whi^ 

35 This smgle equation is commonly refiared to as the IMPES pressure equation. TTie IMPES pressure equation 
may be solved by calling a linear equation solver. Hie pressure solution is then substituted inl» the original set 
of linear equations, and die generalized saturations are con^ted ejqilicitiy. 
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Both the fiiUy in^Ucit method and the IMPES method aim at generating values for the base pressure 
and die generalized saturations at the new dme index n+1 for each cell in the reservoir. However, because the 

IMPES mcfliod is less stable than the My in^licit mctfiod (FIM), die timestep ^^^^ used in die IMPES 

mediod is generally significantly smaller dian die timestep used in die fully iixq>licit mediod. While die 

CE 

5 single-timestep computational effort of IMPES is much smaller dian die single-timestep computa- 



tional effort ^^^^ of die fully implicit mediod, it is quite often die case diat die ratio of timestep 

CE 

sizes is larger than die ratio ^^^^ of con^utational efforts. Thus, any advantage gained by die single- 
timestep efficiency of die IMPES mediod is counteracted by the necessity of performing a large number of IM- 
PES timesteps to cover a timestep of die fully implicit method. 

10 The IMPES mediod is one medxod in a general class of mediods commonly referred to as sequential 

methods. A sequential mediod involves a two-step procedure: a first step in which unknown pressures are de- 
termined, and a second step in which complementary unknowns (i.e. unknowns other than pressure) are deter- 
mined using the pressure solution obtained in the first step. 

Anodier sequential method, commonly referred to as die total velocity sequential semi-inqplicit 

15 (TVSSI) method has received significant use since it was originally developed by Spillette et al. circa 1970. 
The TVSSI mediod is described in die following paper by Spillette, A. G., Hillestad, J. G., and Stone, H. L.: "A 
High-Stability Sequential Solution Approach to Reservoir Simulation," SPE 4542 presented at die 1973 SPE 
Annual Meeting, Las Vegas, Sept. 30-Oct 3. This paper is hereby incorporated by reference. 

Similar to die IMPES method, the TVSSI method has die advantage of reduced con^tational effort 

20 per timestep as compared to die fidly impUcit mediod. However, die TVSSI mediod is far more stable dian die 

IMPES mediod. The increased stability inq)lies diat die timestep ^^rssi of die TVSSI mediod may be signifi- 

candy larger dian die IMPES timestep , The TVSSI mediod comprises two major steps: (i) solving die 

IMPES pressure system; and (ii) solving a set of coiqiled saturation equations for die generalized saturations. 
Since the IMPES pressure equation involves a single unknown (i.e. pressure) at each cell, step (i) requires sig- 
25 niftcandy less work than solving die set of fully implicit equations. In addition, since the set of coupled satura- 
tion equations does not have the elliptic nature of the IMPES pressure equation or die set of fidly implicit equa- 

dons, die saturation sohition converges rapid^. Overall* die single-timestep computational effort for 
die TVSSI method is typically a half to a fifdi diat of die fully inq)licit confutations. 

The TVSSI mediod is not as stable as die folly implicit mediod. In some problems, the ratio ^^^ssr 



CE 

30 of timestep sizes is larger dian die ratio of conqputational efiforts. In odier words, die single timestep 

computational efficiency of die TVSSI mediod relative to die fidly implicit mediod is more dian offset by die 



\ 
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necessity of perfonning multiple timesteps of the TVSSI method to cover a timestep of the fully implicit 
method. 

Overall, the My implicit method seems to be more desirable than the TVSSI method, in part because 
it is more trouble-free. However, the total velocity equations contain a certain powa that enables the success, 
albeit not universal, of the TVSSI method. This power has yet to be Mly appreciated and harnessed. Tims, 
there exists a need for a reservoir simubtion method which may more effectively capture this power inherent in 

the total velocity equatiom. 

One prior-art mdhod used to lower the cost of reservoir simulations is the so called adaptive impKcit 
method (AIM). The adaptive implicit method is based on the recognition that the impUcit fomiulation is re- 
quired at only a fraction of the cells in the reservoir model. If the impUcit formulation can be applied only 
where it is needed, with the IMPES formulation being used at the remaining cells, significant reductions in 
computational effort may be obtained. Ibe adaptive implicit method deteimines dynamically which cells re- 
quire implicit formulation As the simulation progresses in time, a particular cell may switch back and forth 
between IMPES finmulatian and implicit formulatiMi. 

In a rebted prior-ait method, referred to as static variable impUcitness. the assignment of IMPES or 
implicit formulation to each ceU in the reservoir remains fixed through the simulation. 

Although the adaptive impUcit method and variable implicit method are computationaUy more efficient 
than the My implicit mediod, they are still significantly time consuming. Ilus, there exists a need for im- 
proved metfiods for performing adaptive and variable in^Iicit reservoir simulations. 

SUMMARY OF THE INVENTION 

nie present invention comprises a method for performing reservoir simulation by solving a mixed im- 
pUcit-IMPES matrix (MIIM) equation. The MIIM equation arises from a Newton iteration of a variable impUcit 
reservoir model TTie variable implicit reservoir model comprises a plurality of cells including both implicit 
cells and IMPES cells. Tlic MIIM equation mcludes a scalar IMPES equation for each of the IMPES cells and a 
set of inpUcit equations for each of the fat^licit cells. 

In one embodiment, the method for performing reservoir simulation comprises: (a) comtmcting a 
global IMPES pressure matrix equation from the MIIM equation; (b) determining coefficients for a set of satu- 
ration equations at the implicit ceUs by using a total velocity constraim at the impUcit cells; (c) solving the 
global IMPES pressure matrix equation for pressure changes; (d) computing first residuab at the impUcit ceUs in 
response to the pressure changes; (e) solving the set of saturation equations (f«med from the coefficients and 
first residuals) for samration changes at die impUcit ceUs; (f) computing second residuals at the impUcit cells 
and at a subset of tiie IMPES cells fliat are in flow commmiication with any of the impUcit ceUs in response to 
the saturation changes. Steps (b) timmgh (f) may be repeated mitil the second residuals satisfy a convergence 
condition. A final sohition estimate may be computed for the MUM equation from tiie pressures changes and 
the saturation changes after the convergence condition is satisfied. Tbc final sohition estimate may beusedby a 
reservoir simulator to determine behavior of the reservoir model at a fadne discrete time vahie. 

Tie global IMPES pressure matrix equation may be amstmcted firan the MIIM equation by (i) ma- 
nipukting the set of impUcit ceUs at each impUcit ceU to generate a corresponding IMPES pressure equation. 
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and (ii) concatenatiDg the IMP£S pressure equations for the IMPES cells and the IMPES pressure equations for 
ibe iiiq)licit cells. Note the IMPES pressure equations for liie IMPES cells are provided by the MUM equations. 

In a second embodiment, die method for perfoxming reservoir simulation conqnises: (a) constructing a 
5 global IMPES pressure equation from the MUM equation; (b) solving the global IMPES pressure equation for 
pressure changes; (c) con^ting first residuals at the rn^licit cells in response to the pressure changes; (d) de- 
termining improved saturations and improved pressures by perfonning one or more iterations with a selected 
preconditioner at the implicit cells; and (e) computing second residuals at the implicit cells and at a subset of the 
IMPES cells that are in flow communication with any of the implicit cells in response to the improved satura- 
10 tions and inq>roved pressures. Steps (b) through (e) may be repeated until a convergence condition based on the 
second residuals is satisfied. A final solution estimate for the MUM equation may be computed fiom the pres- 
sure changes, inqnoved saturations and improved pressures after the convergence condition is satisfied. The 
final solution estimate may be used to determine behavior of the reservoir model at a future discrete time value. 

15 In a tiurd embodiment, the method for performing reservoir simulation con^rises: (a) constmcting a 

global IMPES pressure equation from the MIIM equation; (b) solving the global IMPES pressure equation for 
pressure changes; (c) conoputing first residuals at the implicit cells in response to the pressure changes; (d) 
solving an iiiq>licit system comprising tiie set of urqslicit equations associated with each of tiie implicit cells for 
improved saturations and improved pressures at the implicit cells using the first residuals at the in^)licit cells; 

20 and (e) computing second residuals for a subset of the IMPES cells which are in flow communication with any 
of the implicit cells. Steps (b) through (e) may be iterated until a convergence condition is satisfied based on the 
second residuals. The final solution estimate for the MIIM equation may be computed based on the improved 
saturations and in^roved pressures after tiie convergence condition is satisfied. In solving the implicit system, 
cell pressures for fringe IMPES cells (i.e. the IMPES cells which are in flow communication with any implicit 

25 cell) are held fixed at those values determined in the pressure solution of step (b). 
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BRIEF DESCRIFnON OF THE DRAWINGS 

Abetterundeistandmgofteprcscminveationcanbeobtamedw 

oflteprefenedembodin^ntsiscon^ 

Fig 1 fllustratcs thestiucture of an inipUcitmalrix «,uati<muscdmicservoksi^^ 
Figs 2A&2Biaustntf»«meenaKKlin««ofalinc«solv« 

Fig. 3 iUustiates a reservoir simulatioa method wMch iovokes a linear solver according to Ae Fesent 

invention; , •f^- • 

Fig. 4 illustrates a reservoir sinndation method which uses total velocity sequential preconditionmg 

according to the present invention; 

Fig 5 iUustrates a partitioning of cells in a variable impUcit reservoir simulation; 

Fig. 6A aiustrates a fust iterative method for solvit* a mixed implidt-IMPES matrix equation accord- 

ing to the present invention; 

Fig. 6B fltastrates a second iterative method for sohnng a mixed impUcit-IMPES matrix equation ac- 

cording to the present invention; 

Fig.7illustratesathird iterative method for solvingamixedimpUcit-IMPES matrix ^^^^^ 

to the present invention. 

While the invention is susceptible to various modifications and alternative forms, specific embodiments 
thereof areshownby way of examplemftedmwings and v^her^nbe described 
stooi however. duutf« drawings and detailed description thereto are not intended to 
particular forms disclosed, but on the contrary, the intention is to cover all modifications, equivalents and alter- 
natives fidling within the spirit and scope of fl« present invention as defined by the 



DETAILED DESCRIPTION OF TOE FREFERRED EMBODIhlENTS 

1. An Implicit Linear Equation SoWer 

He present invention comprises a method for solving an impUcit linear equaticn ^ = C «hich arise, 
fromaNewton iteration of the My impUcit equations. Equation (BIT) above is an example of «^ 
earequation. Matrix A «Ki vector C ate gh^cn. and vector x is to be detetmined. The vector unknown x has the 



0 



!ilpis"av«tor of ceUprtssuresCone pressure per celOandSisavector of cen«^ 
pcrceUforsimuktionswidrMconservedspecies). Given a current estimate L^'-i for the solution of Ae 
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implicit linear equation ^ ^ , the linear solver method of the present invention may be described as fol- 
lows: 

(A) Confute an updated pressure vector ^ using an IMPES pressure equation which is derived 
from the inq)Iicit matrix equation; 

(B) Solve for an updated saturation vector ^ in equations developed using a total velocity conser- 
vation principle; and 

(C) Sappfy the vector - comprising the IKfPES pressure ^"^^ and the updated saturation 
^ to a solution accelerator such as ORTHOMIN or GMRES. 



The updated solution estimate 



- returned by the accelerator forms the basis for the next 



iteration of steps (A) through (C). Steps (A) through (C) are repeated until convergence is attained. 

Let L J represent the intermediate solution estimate after the IMPES pressure vecUv 



is computed. Define vector unknown 



which includes unknown saturation vector 



(It is noted that the pressure vector ^ will not be computed, but its presence here assists in formu- 



lation of the requisite equations). Observe ttut 



(1.1) 



The pressure change P"*^ P"*^ and the saturation change ^ ^ — <S induces a changes in the flow be- 

tween cells. The total velocity equations are used to eliminate the effect of die pressure change P ^P 

on the change in flow. The resulting equations may be solved for the updated saturation ^ 

The linear solver method of the present invention is similar to the combinative mediod in that it in- 
volves a strategy of solving for pressure tet and then for variables other than pressure. 

Each outer iteration of die linear solver method is relatively inexpensive, and success of the method 
hinges on how many outer iterations are needed. The linear solver method is particularly well suited for use 
with the adaptive implicit method (AIM), since the natural way to perform AIM is to begin by solving the global 
set of IMPES equations. 



11 



PCTAJS99/28137 

WO 00/32905 

1.1 Some Theoretical Observations , , 

The linear solver method of the present invention exploits beneficial properties of the total velocity 
equations within a linear equation solver. The linear equation «>lvcrinay be used to soWe an inipUcU hnear 
equatbn ^ = C (When Newton's method is appUed to the My impUcit equations, a whole series of sack 
5 cquationsisBenerated,oneequationperNewtonitera.ioa)TheMowingt^^ 

tivationforthelinearsolvermed^dacomlingtothepresentinvention. The flow velocity between two cells 
is given by the e;q>ression 

v,=A,AO^^ (1.1.1) 

where index denotes a particular phase such as oil. water or gas. Xv is the transmi^ibility-mobility product 

10 forphase^and^^'isthepo.entialdiiferenceforphaseV'be^veenthetwocen^ Ut the subscript b indi- 
cate the b«« phase, ix. the phase whose pressure is solved for in the IMPES pressure equation. Eq. ^ 

be rewritten in a form containing two spatial differences - one that depends 
dependsoncapiUaryprcssure.ie.thedifferenceinpressurebetweenphase V and the base phase b: 

V, = A^,AOj + K^O, - <I>») (1.1.2) 
15 SummingthephasevelocitiesoveraUphasesgivesanexpressionfortotdvelocity^rasfoll^ 

The subscript T denotes a quantity that is summed over all phases V . It can be shown that contimnty con- 
saints force the total velocity to vary substantially less than individual phase velocities. In the extreme case of 
one-dimensional incompressible flow, the total velocity does not vary at aU spatially. 
20 Bysolvingfor jnEq. (1.1 J), and substiuiting the resultant expression into Eq. (l.U). the flow 

velocity may be expressed as 

In anticipation of an iterative method. Eq. (1.1.1) is rewritten in a linearized form: 

25 wherethes«perscriptkistheiterationnmnbera«d&*denotesthecha^^ 
andk+1. Shnilarly.Eq. (1.1.4) may be linearized as 

where 
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** ,and (1.1.7) 

. (1.1.8) 

Finally, an updated phase velocity may be defined as 



(1.1.10) 



5 where 

Eqs. (1.1.5) and (1.1.6) are exactly equivalent If Eqs. (1.1.7) and (1.1.8) are subsdmted into £q. (M.6)» and all 
possible cancellations are performed, £q. (1.1.6) reduces to Eq. (1.1.5). It is noted that Eq. (1.1.6) includes only 

one tenn, i.e. ^ , which depends on the pressure solution ^ . As a result, given an estimate for ^ , Eqs. 

10 (1 .1.6) and (1.1.9) form the basis of a set of equations which invohre only saturation variables ^ . It is noted 

X P =<l)*-<cl>i 

that the transmissibility-mobility products ^ and the capillary pressures ^ ^ are functions of the 

saturation variables. 

\2 Generating Total Velocity Sequential Equations 
IS This section describes the computational steps in generating the total velocity sequential equations 

from the implicit matrix equation where A is a given matrix, C is a given vector, and x is a vector 

unknown comprising cell pressures (one i»ressure per cell) and generalized saturations (M-1 generalized satura- 
tions per cell in a reservoir model with M conserved species). 

Figure 1 illustrates the structure of the in:9licit matrix equation for a reservoir with three cells. How- 
20 ever, the following discussion generalizes to any number N of cells. The matrix A on the left-hand side of the 
implicit matrix equation is an array of submatrices (also referred to herein as blocks) with N block-rows and 2N 

block-cohmms. Eadiof the submatrices ofmatnx Abas M rows and one colunm, where Mis tiie number 

of conserved species. Each of the submatrices ^ of matrix A has M rows and M-1 columns. Thus, nuUrix A 
has NM rows and MM columns. 

P S 
25 The vector unknown x comprises scalar pressures ' and generalized saturation subvectors ' . Hie 

p 

scalar pressure ' is the base pressure at cell i, i.e. die pressure of a predetermined phase at cell L The general- 
ized saturation subvector comprises a set of generalized saturation variables at cell i Therefore, 
vector unknown X has dimension , Vector c, on die right-hand side of die matrix equation, comprisesN 
C C 

subvecton ' . Each subvector ' comprises M known constants. Thus vector C has dimension MN. 
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Each ceU of the reservoir contributes M scalar equations to the nutrix equation. Each block-row of the 
matrix equation sununarizcs the M scalar equations which are contributed by a corresponding cell For exam- 
ple, flie ith block row of the matrix equation, i.e. 

5 suinnuuizes the M scalar equations which are contributed by cell i Equation (U.0) may be equivalently ex- 
pressed in the foim 

which distinguishes (a) the summation term j=i which involves the pressure ' and saturation vector ' for cell 
. i, and (b) the remaining sumiiiation terms which invoWe pressures Itis 

10 noted that the submatrices ^'«' and ^« willbe zero except for those cells j which arc in contact with ccU i 

Each diagonal pressure submatrix may be expressed as the smn of a pressure capacitance subma- 

trix ™ and a pressure (low submatrix 

■ Similarly, each diagonal saturation submatrix ^ may be expressed as the sum of a saturation capacitance 
15 submatrix and a saturation flow submatrix : 

OffMiiagonal pressure submatrices and saturation submatrices ^« relate entirely to flow. each off- 

diagonal pressure submatrix ^'^ may be equated to a conesponding pressure flow submatrix . and each 

off-diagonal saturation submalrix ^ may be equated to a corresponding saturation flow submatrix ^ . 
20 Equation (1.2.1) may be rewritten in a form which distinguishes between capacitance and flow contri- 

buticMis: 



(1.2.2) 



25 



The flow submatrices obey die foUowing relations: 

M , (1.23) 
yw . (1^.4) 
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5 



15 



Thus, the pressure flow submatrix ^ amy be computed by adding the off-diagonal pressure submatrices ^' 
in die ith block-column of matrix A, and negating the resultant sum. Similarly, the saturation flow submatrix 

^ may be conqnited by adding the off-diagonal saturation submatrices ^' in die ith block-column of matrix 
A, and negating the resultant simt 



The Volume Balance Equation 

The volume balance equation combines the M scalar equations at each cell into a single scalar equation 
in such a way that die saturation capacitance disappears. This is accon^lished by detenoining multipliers as 
10 follows. The first step in die detennination of multipliers is to determine die saturation cs^itance coeflicients 
according to die relation 



(L2.5) 



An A^^l vector^' is determmed by solving the linear system given by 

CsuM^^O^ (1.2.6a) 

(1.2.6b) 

where the superscript T denotes ^e matrix transpose operation, and ^ is a vector consisting entirely of ones. 

The con^onents of vector are the multipliers which are used to combine the M scalar equations at cell L 
Since Equation (1.2.6a) conq>rises M-1 scalar equations in M unknowns, an additional constraint is needed to 
obtain unique solutions for die multipliers. Eq. (1.2.6b) is one possibility among many. Another possibihty is 
20 to specify one of die multqsliers, reducing die number of unknowns by one and diereby reducing the computa- 
tional requirement 

The volume balance equation is obtained by pre-multiplying Eq. (1.2.1) by ' . The resulting equa^ 



tionis 



0-2.7) 



25 where 



Bsa=M J Asa = M jF^a ^ ^^2S) 

Bpg^MjApg^ (1^.10) 



IS 
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Bgg-MjA^^ (1.2.11) 

•Die IMPES pressure equation may be obtained ftom Equation {12.1) by evaluating pitssures at intennediatc 
iteration level ^"^^ and saturations at the old iteration level n. Thus, the IMPES pressure equation is as 



follows: 



M >' . (1.2.13) 



The Velocity Equation 

In one prior ait method, ie. the total velocity sequential method, pressures and saturations are com- 
puted according to the following strategy: (a) pressures are computed using Eq. (1.2.13); (b) total velocities are 
computed based on these pressures; and (c) saturations are computed while holding fixed the total velocities. 

Since velocity relates to flow between connections, rather than conservation at a cell, the equations of 
the present invention require a different constraction. Let be the vector of flows from cell i to cell j defmed 
by 

where the M coir^onents of '^^^ represent the flows of each of the conserved species from cell i to ceU j; ^« 

and arc vectors; and and are ^^(^"l) matrices. These terms can be extracted from 

the original matrix equation, as foUows. The two off-diagonal (i.e., j) terms are 

Fpj = ^P\} ^ (1.2.15) 

(12.16) 

Hie corresponding flows fixan cell j to cell i obey Ac relation 

p^^^ffi (1.2.17) 
As a result, the diagonal (Le., i) terms can be obtained from Ae equations at die connected cells, leading to 

Ppi-^^iyi ^ (12.18) 

(1.2.19) 

The view from cell i of the total volumetric flow from cell i to cell j is given by 

TTie view from cell j of the total flow from ceU j to ceU i, in addition to having an opposite sign, has a different 
magnttiiHe becausc its vcctor of multipliers is different, i.e. 
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'"J^ . (1.221) 
The total flow, as viewed fiom cell i, is given by the following oqjiession, which is obtained by multiplying 

A/' 

(U.14)by . 

F» = + + F^Pj + F«,S, + F^Sj 



(1222) 



where 



F^''=MjF«^ 

» 

= mJf^ 

I 

By solving the IMPES pressure equation (1.2.13), pressures ' at the intermediate iteration level 

^'^y^ are obtained These pressures may be substituted into Equation (1.2.22) giving an expression for the 
total velocity after the IMPES pressure solution: 

F/ +F^r^ ^F^P;*>i^F^S;^F^S'; ^j^23) 

In the discussion to follow, a set of equations will be developed which enable the computation of a new set of 

saturations ' . Let ' be the pressures which correspond to die new saturations and are used to compute 
flow from cell i to cell j. (lliese pressures will not be computed, however it is convenient to include diem here 
to aid the development of the requisite equations.) Ibe total velocity corresponding to die new pressures and 
saturations is given by 

F/ = Fr + F^If + F^Pf + F^sr^ + F^sr^ 

Thus, cell i's view of the change in total velocity is obtained by subtracting Equation (\22^) finm Equation 
(1^.23): 

+F^(sr^-s:)+F^(sy^-s^) 



(122S) 



gp9 

This total velocity chaqge ^ is set equal to zero, yielding 

F?M-pr''VF^{p!-pr'') 

+F^(sr^ -5;)+F4(s7'i -5;)= o 



F" 

Note that tiiepiessuie at cell i in (12^6), ' ,niay vaiy withj. 



(1226) 
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Writing the balance equation (1.2.2) in terms of the desired iteration levels and rearranging yields 

1 *>i J J^^ l^J J {1227) 

Equation (1.2.26) may be used to eliminate the pressure difference ^ ^ &om Equation (1.2.27). How- 
5 ever, it would be advantageous if Equation (1.2.26) could be used to eliminate the pressure difference 

^ - if Equation (1227) at Ae same time. The most likely conditions that would permit the elimi- 

nation of both pressure differences is to have 

^H-'^fi^ {122%) 

^TPi- ^7r,\ (1.2.29) 

10 These relatiom arc approximately trae, and it win be assumed that if die pressure ^ is 

eliminated the other pressure difference ^ " ^ wttl disappear as well. In effect, when equation (1.2.22) 

IS used to eliminate ^^^^^ it is assumed that this elimination simultaneously eliminates this con- 

F Ip^^P""^] 

nection'scontributiontotheproduct'^'"^ ' ' The resultant equation is 
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1.3 Solution oftheInq>licit Matrix Equation 

The equations developed above arc used to constract a linear equation solver. It is assumed that a res- 

ervoir simulator executing a fuUy inqslicit simulation generates an implicit linear equation of the form 

Hie reservoir simulator provides the matrix A and vector b as input data to the linear equation solver of the pre- 

20 sent invention. Tlie linear equation solver returns an estimate for the solution = ^ * to the reservoir simu- 
lator. TTie linear equation solver enyloys an iterative metbod according to the present invention for solving the 

implicit linear equation. Each itwation operates on a current estimate and generates an updated estimate 

The sequence of estimates , ^\ ... generated by the linear equation solver converge to 
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the solution ^ of the implicit linear equation. Given a current estimate L -I, where ^ isavector 

of cell piessures (one pressure per cellX and is a vector of saturations (M-1 saturations per cell for reservoir 
models with M conserved species), a generic iteration of the linear solver method includes Ae following steps. 

5 1. Construct the IMPES pressure equation (1^.13). This is achieved by performing the 

column summation indicated by Eq. (I JS.5}. In other words, for each diagonal satu- 

A A 
mtion subroatrix ^ of matiix A, the diagonal saturation submatrix is added to 

A 

the off-diagonal saturation submatnces in the same block-column, thereby gen- 

C 

erating a corresponding saturation capacitance submatrix . The saturation ca- 

10 pacitance submatrix relates to the accumulation of each of the species in cell i. 

Usmg the capacitance submatiices, the well-known IMPES reduction is performed. 
The IMPES reduction results in matrix equation (1.2.13) which involves only the 
pressure variables. 

15 2. Construct the saturation equation (1.230) as described above\ 

3. If necessary, compute die underlying implicit equation residuals. 

4. Based on the current in^ticit equation residuals, compute the IMPES pressure equa- 
20 tion residuals. 

5. Solve the IMPES pressure equation (1.2.13) for updated pressures ^ ,andcom- 
pute a corresponduig set of pressure changes ' ' , where ' denotes cell 



25 



pressures at the beginning of the current iteration. 
6. Update the in^licit equation residuals for the pressure changes computed in step 5. 



7. Solve ihe saturation equation (1.2.30) for updated saturations ' , and compute a 



30 



corresponding set of saturation changes ' ' . 



8. Update the implicit equation residuals for die saturation changes computed in step 7. 



19 



wo 00/32905 



PCT/US99/28137 



9. Combine Ac pressure and saturation changes of steps 5 and 7 into a con^osite solu- 

tioachange *- 

10. Feed this solution change ^ to a solution accelerator such as ORTHOMIN or 
GMRES to accelerate conveigence. 



11. Update the implicit equation residuals based on the solution estiniate 
returned by the accelerator. 



Steps 4-1 1 are repeated untU convergence is attained. 

Figures 2A & 2B iUustrate one embodiment of the linear solver mediod accoiding to the present inven- 
tion. The linear solver method shown in Figures 2A & 2B may be implemented in software on a computer sys- 
tem. The linear solver method is typically invoked by a reservoir simulator also in^lemented in software. The 
reservoir simubtor provides the linear solver method with an inq)licit matrix equation = * which resuhs 
from a Newton iteration on the fiiUy implicit equations. The linear solver method comprises the following steps. 

In step 110. a global IMPES pressure equation is constructed from the waphcit matrix equation 
Ax = b^ The global IMPES pressure equation may be constructed as described above in the development of 

IMPES pressure equation (1 J1.13). 

In step 120, Uie global IMPES pressure equation is solved to determine an improved estimate of pres- 
sure at a plurality of cells. In the preferred embodiment, the pluraUty of cells include all the cells of the reser- 
voir. In another embodiment, the plurality of cells may represent a subset of the cells of the reservoir. 

In step 130, residuals of the implicit matrix equation are updated based on the improved estimate of 
pressures. 

In step 140, a conqilementary matrix equation is constructed in terms of unknowns otiier than pressure. 
The convlemcntary matrix equation is constructed from tiie irap^cit matrix equation based on the constraint of 
preserving total velocity between cells. For exan^le, the con^lementary matrix equation may be saturation 
equation (1.2.33). 

In step 150. the complementary matrix equation is solved in order to determine an inq>roved estimate 
of the unknowns other than pressure at each cell of the reservoir. 

In step 1 60, the residuals of the inq>licit matrix equation are updated based on the improved estimate of 
the unknowns other than pressure. 

In step 170, a conqwsite sohition change which conqirises a first change in pressure associated witfi Ac 
improved estimate of pressures determined m step 120 and a second change in the unknowns other than pressure 
associated widi die inqnoved estimate of die unknowns other than pressure. 
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The composite solution change is treated as the output of a preconditioner. In step 180» the composite 
sohition change is provided to an accelerator such as, e.g., GMRES or ORTHOMIN, in order to accelerate con- 
vergence of the solution. 

In stq> 190, the solution accelerator generates an accelerated sohition change. 
5 In step 195, the residuals of the in^licit matrix equation are updated based on die accelerated solution 

change. 

In step 200, a test is performed to determine if a convergence criteria has been satisfied. If the conver> 
gence criteria is not satisfied, another iteration of steps 120 through 195 is perfomied. If the convergence crite- 
ria is satisfied, a final solution estimate is conq>uted based on the accelerated solution change and a previous 
10 solution estimate as mdicated by step 202. 

In step 205, the fmal solution estinute is applied to predict the behavior of reservoir fluids at a future 
time value. 

In one embodiment of the linear solver method, the complementary matrix equation is a saturation ma- 
trix equation such as equation (1.2.33), and the unknowns other than pressure are saturations. 
15 In another embodiment, the unknowns other than pressure comprise one or more variables such as, e.g., 

saturation, mole fraction, mass, energy, etc. 



Figure 3 illustrates the structure of a reservoir simulator method which invokes the linear solver 
20 method as described above. In step 310, the reservoir simulator formulates a set of finite difference equations 
which describe a generalized timestep in the time evolution of fluid properties in the cells of a reservoir. In step 
320, the reservoir simulator performs one or more Newton iterations in order to solve the finite difference equa- 
tions for a single timestep. The solution of the finite difference equations defines a pressure and one or more 
complementaxy unknowns for each cell in the reservoir at die next discrete time level 

25 

Each Newton iteration con^srises the following steps. In step 320A, a linear approximation is con- 
stmcted for each of the non-linear terms in the finite difference equations. In step 320B, an implicit matrix 
equation is constmcted based on the finite difference equations and the linear approximations. In step 320C, the 
implicit matrix equation is solved using the linear equation solver method discussed above in connection with 
30 Figures 2A&2B. 

By performing a series of timesteps as described above, the reservoir simulator may predict the behav- 
ior of the reservoir fhiids. 

35 1.4 A Preconditioner for Solving die Implicit Matrix Equation 

The present invention also con^mses a preconditioning method for solving the implicit matrix equation 

^ ^ . The preconditioning method has performed effectively in a variety of problems. 
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Given a cuirent estiniate J for the solution to the in^Ucit matrix equation, is a 

vector of cell pressures and is a vector of cell satuiations. the preconditioning method of the present inven- 
tion con^rises Ifae following steps: 

(1) Solve the IMPES pressure equation for an updated pressure vector and compute pressure 
change vector ^ ^ ; 

(2) Update the implicit equation residuals for the pressure change ^ ^ ; 

(3) Solve saturation equations (1^33) for updated saturation vector ^"^^ , and compute saturation change 
vector and 



(4) Update the implicit equation residuals for the saturation change ^ 



Ite composite solution change L^" J is supplied to a solution accelrator such as 

Ordiomin or GMRES. 

Any suitable method can be used to solve the IMPES pressure equation. The saturation equations tend 
to be easy to solve, in the sense that an iterative solution of saturation equations converges rapidly. This sug- 
gests use of a simple prcconditioner such as diagonal scaling or ILU(0). ILU(0) was used in the tests described 
below. 

TTie preconditiMing method of fte present invemicm differs from 4e Constrained Pressure Residual 
Method (Wallis, J. R., Kendall, R. P., and Uttle, T. E.: "Constrained Residual Acceleration of Conjugate Re- 
sidual Methods," SPE 13536 presented at the SPE 1985 Reservoir Simulation Symposium, DaUas. Texas, Feb- 
ruary 10-13, 1985) in at least two ways. First, the preconditioning method of the present invention obtains the 
pressure equation using the true IMPES reduction. Wallis et aL perform a reduction direcUy on the implicit 
equations. Second, the prcconditioner method of the present invention solves the total-velocity saturation equa- 
tions. Wallis etaLperfoim a single itention on the implicit equations using a precondition 

system ILU(0). 

Test Restths and Discussion 

Tbe prcconditioner method has been tested on a handful of matrix equations. Table 1 below summa- 
rizes the results. The convergence criterion used was a 0.005 reduction in the residual L2 nomt Case 1 was 8 
variant of the first SPE comparison problem (Odeh, A. S.: Comparisons of Solutions to a Three-Dimensional 
Black-Oil Reservoir Siimilation Problem," JPT 33, January 1981. 13-25), with the wells being treated as flowing 
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against constant pressure. Case 2 was the first Newton iteration of the first timestep of the ninth SPE compari- 
son problem (Killough, J. E.: '*Ninth SPE Con^arative Solution Project: A Reexamination of Black-Oil Simu- 
lation,** SPE 291 10 presented at the 13tfa SPE Symposium on Reservoir Simulation, San Antonio, Texas, Febxu- 
ary 12- 1 6 J 995), with a one day timestep. Case 3 was die same as case 2, with the timestep size increased to 50 
days to make tiie problem more diflBcult. Case 4 was the same as case 3, but for the second Newton iteration. 
Case 5 was from a 2400-ceU, two-hydrocaxbon CQnq}onent, steam injection model. Case 6 was fit>m a 5046- 
cell, seven-hydrocarbon conq>onent plus water conq)ositional model 

The logical conq)arison to make is to the Consttained Pressure Residual (CPR) Mediod. In these problems, 
the new method took either the same number as or somewhat fewer outer iterations than CPR It required on 
average a little less than two saturation iterations' per outer iteiadon. The resulting conqnitational woik required 
was probably somswbat less than dial requured by CPR's single reduced-system ILU(0) iteration. 



Table 1 — Performance of Two-Stage Preconditioner 
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1 


2 


4 


3 


8 
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Total Velocity Preconditioning in a Reservoir Simulator 

Figure 4 illustrates a reservoir simulation method which uses total velocity sequential preconditioning 
according to the present invention. The reservoir simulation method coxrqnises the followmg steps. In step 410, 
the reservoir simulator formulates a set of finite difference equations which describe a generalized timestep in 
the time evolution of fluid properties such as pressure, saturation, etc. for each cell in the reservoir. 

In step 420, the reservoir simulator solves die finite difference equations by performing one or more 
Newton iterations. The solution of die fmite difference equations specify the value of pressure and complemen- 
tary unknowns (i.e. unknowns other tiian pressure) for each cell at the next time level For each Newton itera- 
tion, die reservoir simulator 

(a) Constracts a linear approximation for each of the non-linear terms in the finite difference 
equations as indicated by step 420A; 

(b) Constructs an inq)licit matrix equation based on the finite difference equations and die 
hnear approximations as indicated by step 420B; and 
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(c) solvesthein>plfcitinatrface«p»tionby(cl)constmc^^ 

to teons of imknowiu other than pn«sure. and (c2) solving the complementary matnx 
«p««ion for the mtoovvm other than pr«»«rc as indicated by step The comply 
^ matrix equation is constructed using a constraint of conserving total veloaty 
between cells. 

By performing a succession of timcstcps. Le. by repeatedly solving the finite difference equations, the 
.i^evolullTofprcssureandthccomplen^unknow^ TTus information may be used. 

c.g.. to guide toe development andmanagementofaphysical reservoir suchasanoa field. 

LSLinearSolvetsforVariablelmplicitandAdaptivelmplicitSimulations 

ITU, present invention comprisesamcthod for so^ingthcmatxixcquadonsw 
pHcit^d adaptive implicit reservoir simulations. As discussed above, the fi% in^Ucit formulat.onre^^ 
^^Itly more computational effort per timestep than the IMPES fonnulation. However, the larg. 
Zlsltmaybeusedw«h*eMyh.^Ucitfen».a.onoftenmare^o£.ets^addrn^ 

'""^ '""^e nonlinearity of the My implicit formuUtion requires an i^^^^ 
EachNewtoniterationgencratesamatri.equationrefe.redtohereinas.hei^^^^ 

timestep of the My impHcit formulation requires the solution 

Dlains the large computational effortof the folly implicit f^^ . • • , 

' one priZlthod used to lower the cost of reservoir simulations is the so called adapnve «nph«t 
„^od (AIM) m adaptive implicit med.od is based on the recognition .hat the implicit fommlation « re- 
X^a^^a .acJof the celU^the reserve, model. If ^e^UcUform^^ 

it is uLed, With the IMPES formuUtion being used at the remains 
:^u..oaale«brtma,beob.a^ ^ ...p^ in^Hcit me^ determines dynamicaUyw^^^^^^^ 
JTimplicitformulation. M the simulation progre«es in time, a p«dc«h. cell may swrtch b.^ 
between IhD»ESfotmuhition and implicit formulation. ^.^fxMP^ot 
m a «Uted prior-art method, referred to as static variable impUcrmess. the assrgnment of IMPES or 
implicitformulationtoeachceUmthereservoirremainsfixedthroughthesimul^ion. 

Mvariableimplicitandadaptiveimplicitreservoirsimul^^ 
describetheimplicitcellsand*elinearIMPESequa.ion,which^^ 
««con^sitesys^ofe^tion,ftom.n^censisnon.^arand.e<^aNewt^^^^ 
compI^^issoWedinaseriesofNewtoniteratio™,. Each Newton iteration results ma mncednnpUcU- 

IMPES matrix equatioa , 

solution of the mixed implicit-IMPES matrix equation poses a challenge to a linear equaUon solver 
TOssectiondescribes two related methods according to theprcsentinven*^ 

solving the mixed implicit-IMPES matrix equation. 

Whenvariableimplicitnessisusedinareservoirsimulation.only^ 
pe^ of the cells are treated implicitly. As shown in Figure 5. the impUcit cells tend to appear as smaU ^ 
L(lg.ishu«lsA,B.C«u,D)inamuch larger IMPES ocean E. At the IMPES ce.^^^^ 
^ 1 be solved fbr. and correspondingly there is a single equation .o be solved. At the unphcrt cells, the 
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number of unknowns is equal to the number of components (such as, e.g., oil, water and gas) being used in the 
model. 

Solving die Mixed Inq>licit-IMPES Matrix Equation: Method 1 
5 This section presents a fust Imear solver mediod according to the present invention for solving the 

mixed implicit-IMPES matrix equation = C , "Hic vector unlmown x comprises a set of cell pressures ' 

c 

(one pressure per cell) and a set of saturations ^ (M-1 saturations per cell for simulations with M conserved 
species). The matrix A and the vector C are supplied to the linear solver method as inputs by a reservoir simu- 
lator. The linear solver method generates an estimate for the solution ^ ^ to the mixed inqjlicit-IMPES 
10 equation. The linear solver metliod comprises an iterative procedure. Each iteration of the linear solver method 

operates on a cuirent solution estimate ^ and generates an updated solution estimate ^ . The sequence of 

solution estimates , ^\ -^^j ... converges to the solution ^ of the mixed implicit-IMPES equa- 

tion. The linear solver method en^)loys a convergence criteria to determine when iterations should temunate. 
Each iteration of the linear solver method comprises the following steps. 

15 

1. Construct a global IMPES pressure matrix equation from the mixed implicit-IMPES ma- 
trix equation. The global IMPES pressure matrix equation con^)rises one scalar IMPES 
equation per cell of the reservoir. The nuxed implicit-IMPES equation already specifies 
the scalar IMPES pressure equation for each of die IMPES cells. At each of the implicit 
20 cells, a scalar IMPES pressure equation may be generated by combining the in^hcit 

equations according to the procedure described above in die sections entitied "Generating 
Total Velocity Sequential Equations'* and The Volume Balance Equation". 



25 



2. Compute the coefficients for die saturation equations (1.2.30) at the in^licit cells. 

3. Solve tile global IMPES pressure matrix equation for mtezmediate pressures ^ ,Le. a 
single intennediate pressure at each cell in the reservoir, and confute pressure changes 

^ " ' based on die intermediate pressures ' and pressures ^ prevailing at 
the beginning of the iteration. 

4. Update implicit equation residuals at die implicit cells based on the pressures changes 

i i computed in step 3. 

5. At the implicit cells, solve for inp:oved saturations ' in saturation equations 
35 (1.2J0) which are derived using a constraint oftotal velocity conservation between cells. 
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6 Update implicit equation residuak at the impUcit cells and at the fringe of IMPES cells 
that arc in flow communication with the impUcit cells based on the satttration solutions 
obtained in Step 5. 



7, Detennine if a conveigence condition is satisfied. 

steps 2-6 axe repeated until the «mverBence condition is satisfied. Note that at the end of step 6. the 
only cells where the residuals faU to meet the convergence criteria are the implicit cells and the fiinge of IMPES 
cells m flow commnnicadon with any hnplicit ceU. Tl»e residuaU at the IMPES cells outside the fringe sttU are 
at the values they had foUowing the IMPES solution. Tte means *at ORTHOMIN or GMRES compulsions 
need be applied only at Aese cells, Le. at the implicit cells and fringe IMPES cells. 

Fig«e 6A illustrates the first method for solving the mixed impUdt-IMPES matrix equation accordmg 
ro the present mvemion. n« mixed impUcit-IMPES matrix equation specifies a set of implicit equations for 
each implicit ceU aada single scalar IMPES pressure equation for each IMPES cell. 

ta step 1010. a scalar IMPES pressure equation is constmcted for each of the impUcit cells. The scalar 
IMPES pressure equation for an implicit ceU is generatedby forming a linear combination of the implicit equa- 
tions which correspond to the implicit cell 

In step 1020. a global IMPES pressure matrix equation is constructed by concatenating the scalar IM- 
PES pressure eqwtions for the implicit cells with the scalar IMPES pressure eq^^^^ 

scalar IMPES pressure equations for the IMPES cells are provided by the mixed impUcit-IMPES matmc equa- 



tion. 



to step 1025. coefficients for a set of sawration equatiom are determined at the implid^ 

total velocity constraint at the implicit cells. 

fa step 1030. the global IMPES pressure matrix equation is solved for pressure changes. In step 1035. 
tilie residuals at the hnpUdtcelk are computed inresponse to the pressure changes determined in step 

fa step 1040. the set of saturation equations are solved at the impUcit cells. TTie set of saturation equa- 
tions are fonned usmg the coefficients (determined in step 1025) and the residual computed in step 1035. 

fa step 1050. hnplicit equation residuals (i.e. residuals at the impUcit cells ami at the fringe of IMPES 
cells that are in flow commmiication with the impUcit celb)are updated inresponse to the saturationd^ 

fa step 1060. a convergence condition is tested based on the updated residuals. If the convergence 
conditionisnotsatisfied.processh«continue,with«Krtherite^^ If the eom«gence condition 

is satisfied, the method termhutfes ami the find sototion esthnate is provided to the caUm^ 
crally a reservoir sumOator. When the convergence condition is satisfied, it is assumed that the solution to the 
nrixed impttdt-IMPES equations has been detemuned with acceptable accuracy. TTie fmal solution estmiate 
comprises a set of converged saturations and pressures which are used by the reservoir simulalor in modehng 
characteristics of tiie reservoir. 

Solving the Mixed Implicit-IMPES Matrix Equation: Mefliod 2 



26 



wo 00/32905 



PCT/US99/28137 



This section presents a second linear solver method according to the present invention for solving the 

mixed inq>licit<-IMPES matrix equation = ^ . Each iteration of the second linear solver method comprises 
the following steps. 

5 1. Ck)nstxuct a global INfPES pressure matrix equation fom die mixed imp^^^ 

trix equatioa The global IMFES pressure matrix equation comprises one scalar IMPES 
equation per cell of the reservoir. The mixed inq>licit-IMPES equation abready specifies 
the scalar IMPES pressure equation for each of the IMPES cells. At each of the implicit 
cells, a scalar IMPES pressure equation may be generated by combining the unplicit 
10 equations according to the procedure described above in die sections entitled "Generating 

Total Velocity Sequential Equations" and *The Volume Balance Equation". 

2. Solve the global IMPES pressure matrix equation for intemiediate pressures ' , ie. a 
single intermediate pressure at each cell in die reservoir, and compute pressure changes 

IS ' ' based on die intermediate pressures ' and pressures ' prevailing at 

the beginning of the iteration. 

3. Compute implicit equation residuals at the implicit cells based on the pressures changes 



20 



computed in step 2. 

4. At the in^licit cells, solve for unproved saturations ' and second intennediate pres- 
p'**H 

suies ' by performing one or more iterations with a selected preconditioner such as, 
e.g.,ILU(0). 

25 5. Update implicit equation residuals at the implicit cells and at the fringe of IMPES cells 

that are in flow communication with the inq)licit cells based on the inq^roved saturations 
and second intermediate pressures obtained in step 4. 



6. Detemune if a convergence condition is satisfied. 

30 

Steps 2-6 are repeated until the convergence condition is satisfied. Note that at the end of step 5, die 
only cells where the residuals fail to meet die convergence criteria are die in^ilicit cells and die fiinge of IMPES 
cells in flow communication with any implicit cell. The residuals at die IMPES cells outside die fiinge stiU are 
at die vahies they had following the IMPES solution. This means diat ORTHOMIN or GMRES computations 
35 need be applied only at these cells, Le. at die implicit ceils and fiinge IMPES cells. 
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Figure 6B fllustrates Ac second metf^od for solving the nnxed in^plicit-IMPES n^ equation ac- 
cording to the present invention. 'H.e mixed implicit-IMPES mattix equadon specifies a set of i«plK»t equa- 
tionsforeachinipMtceUandasinglescalarlMPESpressureequationfore^AI^^ 

Instep 1060,ascdarIMPESpn»sureeq«ationiscon«ructedfbr«^ 
IMPES pressure equadon for an implicit cell is generated by fonning a linear combination of the imphat equa- 
tionsvi*ichcorreq>ond to Ae implicit cell 

in step 1065. a global IMPES pressure matrix equation is constructed by concatenatmg the scalar IM- 
PES pressure equations for ti.e impUcit ceUs with the scalar IMPES pressure equation 
scalar IMPES pressure equations for the IMPES cells are provided by the mixed implicit-lMPES matrtx equa- 

Instepl070.thegtobdIMPESpressurematrixequationissolvedforpressurechang^ Instepl075. 
ti^residuals at the implicit cells are computedinresponse to thepressurech^^^^ 

to step 1 080, improved samrations ami improved pressures at the impUcit cells may be detemnned by 
pcrfomtingoneormoreiteradonswithaselectedpreconditionersuchasILUCO). 

In step 1090. implicit equation residuals (i.e. residuals at the implicit cells and at the fhnge of IMPES 
cells that are in flow communication with die impUcit cells) are updated in response to the unproved satimitions 
and improved pressures. 

to step 1095. a convergence condition is tested based on the updated residuals. If the com^ergence 
conditionisnat satisfied, processing continues withanotheriterationofst^ 1070. If the convergence condition 
is satisfied, flu: method terminates and the final sohtion estimate is provided to the calling routine whrch . 
erally a reservoir simulator. When the convergence condition is satisfied, it is assumed that the solut^n to the 
n^ed implicit-lMPES equations has been detemnned with acceptable accuracy. The fmal soU«ion estimate 
comprises a set of comrerged saturations and pressures which am used by .he reservoir simuUtor in modehng 
characteristics of the reservoir. 

Solving the Mixed Implicit-lMPES Matrix Equation: Method 3 

to this subsection, a third method according to the present invention for solving ti.e m«ed m-phcrt- 
IMPES matrix equation is presented. The stn^ture of tins second method 
„„tiroddescn1,ed above exceptinsteps4and 5. h» this thirdtnethod,,^^ 

and 511 respectively. 

411 Solveforsatinations andpressures ^'^^ at ti«unplicit cells while holdtog fixed the pre,- 
s„«s in the surrounding ftinge (of IMPES cells) to ti« values^' determined during the IM- 
PES pressure solution. Any method can be used to generate the solutions for satiuations > 

and pnssures . hut it must be able to deal with the unstructured fomi of the impU^^ 
equations. 
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511. Update residuals in the fringe of IMPES cells. Since the implicit equations have been solved, their 
residuak will satisfy ^ convergence criteria. 

After step SII, only die fringe cells will have residuals that fail to meet the convergence criteria. Again, OR- 
5 IHOMIN or GMRES only need be ^lied to these cells. 

Figure 7 illustrates one embodiment of the diird mediod for solving the mixed implicit-IMPES matrix 
equation according to the present inventioiL The mixed implicit-IMPES matrix equation specifies a set of im- 
plicit equations for each tn^licit cell and a single scalar INfPES pressure equation for each IMPES celL The 
10 enibodiment of Figure 7 comprises the Mowing steps. 

In step 1 1 1 0, a scalar IMPES pressure equation is constructed for each of the implicit cells. The scalar 
IMPES pressure equation for an implicit cell is constructed by forming a linear combination of the implicit 
equations which correspond to the implicit cell. 

In step 1120, a global IMPES pressure matrix equation is constmcted by concatenating (a) die scalar 
15 IMPES pressure equations for the implicit cells and (b) the scalar IMPES pressure equations for the IMPES 
cells. The scalar IMPES pressure equations for the IMPES cells are provided directly by die mixed implicit- 
IMPES matrix equation. 

In step 1 130, the global IMPES pressure equation is solved for pressure changes. 
In step 1 135, the residuals at the in^Iicit cells are computed in response to the pressure changes deter- 
20 mined in step 1 130. 

In step 1140, improved saturations and improved pressures at the implicit cells are determined by 
solving the system ofhaplicit equations associated with the implicit cells while holding fixed the pressures in 
the fringe of IMPES cells which are in flow communication with any implicit celL 

In step 1150, the residuals in the fringe of IMPES cells (which are in flow communication witii any 
25 in^licit cell) are updated. 

In step 1160, a convergence condition is tested based on the updated residuals. If the convergence 
condition is not satisfied, the method continues with a next iteration of step 1 130. If the convergence condition 
is satisfied, iteration terminates and the final solution estimate is returned to the calling routine (e.g. a reservoir 
simulator). The converged saturations and pressures making up die final solution estimate are used by the res- 
30 ervoir simulator in modeling characteristics of the reservoir. 

Observations 

Methods 1 and 2 are less expensive than Mefliod 3 per outer iteration. In **easy** problems, only one it- 
eration may be needed, so Methods 1 and 2 would be preferred In "hard" problems. Method 3 requires fewer 
35 outer iterations. As die problem becomes harder, Mediod 3 becomes preferred. Method 3 effectively requires 
an unstructared implicit equation solver. If such a sohrer is not available, Methods 1 and 2 are much easier to 

40 
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WHAT IS CLAIMED: 

I A method for perfoiming reservoir sinmlation by solving a mixed implicit-IMPES matrix 
(Mm eqaatio... wherein the MUM equation arises ftom a Newton iteration of a variable impKcit reservoir 
nuHiel. wherein the variable implicit reservoir model comprises a phrality of cells mcluding unpUdt cells and 
IMPES cells, vAerein the MIIM eq.«UioniDch«le, afm^tscato 

cells and a first set of implicit equations for each of the impUcit cells, the method comprising: 

a) constructing a global IMPES pressure matrix equation ftom the MIIM equation, wherein said con- 
strocting fte global IMPES pressure matrix equation comprises: 

constructmg a second IMPES pressure equation for each of the impUcit cells ftom the first set 

of implicit equations corresponding to flie implicit cell; and 

concatenating the first scalar IMPES pressure equations for the IMPES celb and the 

IMPES pressure equations for the iaaplicit cells; 

b) determining coefficierts for a second set of saturation equations at the impUcit cells by using a total 

15 velocity constraint at the implicit cells; 

c) solving the global IMPES pressure matrix equation for pressure changes; 

d) computing first residuals at die implicit cells in response to the pressure changes; 

e) solving the second set of saturation equations for saturation changes at the implicit cells, wherein the 
second set of saturation equations are formed with the coeffidents and the first residuak 

20 f) computing second residuals at d>e implicit cells and at a subset of the IMPES cells that are in flow 

communication with any of the impUcit cells in response to the saturation changes; 

g) determining if a convergence condition based on tiie second residuals is satisfied; 

h) repeating b) through g) until die convergence conditio is satisfied; 

0 computing a fmal solution estimate for the MIIM equation ftom the pressmes changes and the satu- 
25 ration changes after Ae convergence condition is satisfied; 

j) applying the final sohition estimate to determine behavior of the reservoir model at a future discrete 

time value. 

2. A method for performing reservoir simulation by solving a mixed unplicit-IMPES matrix 
30 (MEM) equation, wherein die MIIM equation arises fiom a Newton itention of a variable inyUcit reservoir 
n»del wherein the variable implidt reservoir model comprises a phnality of cells hud^^ 
IMPES cells, wherem the MIIM equation mchules a fns. scalar IMPES equation for each of the IMPES cells 
and a set of implicit equations for each of the implicit cells, die metiiod comprising: 

a) constructingaglobal IMPES pressure equation from the MIIM equation, wherein said constructi^ 

35 the global IMPES pressure equation conqmses: 

constructing a second scalar IMPES pressure equation for each of the unplicit cells ftom die 

set of implicit equations conesponding to die implicit cell; and 

concatenating the first scahtt IMPES pressure equation for each of die IMPES cells and die 

second scalar IMPES pressure equation for eadi of die in^Udt cells; 
40 b) solving die global IMPES pressure equation for pressure changes; 

c) computing first residuals at die impUcit ceUs in response to die pressure changes; 

30 
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d) detemiining in^roved saturations and improved pressures by performing one or more iterations with 
a selected preconditioner at the implicit cells; 

e) computing second residuals at die in^licit cells and at a subset of the IMPES cells that are in flow 
communicadon with any of die implicit cells in response to die inqiroved saturations and in^roved pressures; 

5 0 determining ifaconvergaicecanditionbased on die second residuals is satisfied; 

g) repeating b) througih f) until the convergence condition is satisfied; 

h) confuting a fmal solution estimate for the MUM equation irom the pressure changes, improved 
saturations and improved pressures after the convergence condition is satisfied; 

i) applying the final solution estimate to determine behavior of the reservoir model at a future discrete 
10 time value. 

3. A method for performing reservoir simulation by solvmg a mixed in^licit-IMPES matrix 
(MIIM) equation, wherein the MIIM equation arises from a Newton iteration of a variable implicit reservoir 
model, wherein the variable implicit reservoir model con^rises a plurality of cells including implicit cells and 

15 ' IMPES cells, wherein the MIIM equation includes a first scalar IMPES equation for each of the IMPES cells 
and a set of implicit equations for each of the implicit cells, the method comprising: 

a) constructing a global IMPES pressure equation fiom the MIIM equation, wherein said constructing 
die global IMPES pressure equation conqnises: 

' constructing a second scalar IMPES pressure equation for each of the implicit cells from the 
20 set of implicit equations corresponding to the implicit cell; and 

concatenating the first scalar IMPES pressure equation for each of the IMPES cells and die 
second scalar IMPES pressure equation for each of die implicit cells; 

b) solving die global IMPES pressure equation for first pressure changes; 

c) concqniting first residuals at the implicit cells in response to die first pressure changes; 

25 d) solving a system comprising the set of implicit equations associated with each of die implicit cells 

for improved saturations and improved pressures at the implicit cells using die first residuals at the implicit 
cells; 

e) computing second residuals for a subset of die IMPES ceUs which are in flow communication with 
any of the implicit cells; 

30 f) determioiiig if a convergence condition is satisfied based on die second residuals; 

g) rq)eatedly performing b) through f) until the convergence condition is satisfied; 

h) confuting a final solution estimate for the MIIM equation using the uiqsroved saturations and imr 
proved pressures after die convergence condition is satisfied; 

i) applying die final solution estimate to detennine behavior of the reservoir model at a future discrete 
35 time vahie. 

4. The mediod of claim 3, wherein said solving die system for in^roved saturations and im- 
proved pressures at the implicit cells comprises holding cell pressures fixed for a subset of the IMPES ceUs 
which are in flow communication widi any of die implicit cells. 

40 
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5. A method for peifomiing reservoir simulation by solving a mixed ixi^licit-IMPES matrix 
(MIIM) equation, wherein the MHM equation arises from a Newton iteration of a variable is^licit reseivoir 
model, wherein the variable implicit reservoir model comprises a plurality of cells including in^licit cells and 
IMPES cells, wherein die MUM equation includes a first scalar IMPES equation for each of the IMPES cells 

5 and a set of implicit equations for each of the implicit cells, the method conqsiising: 

a) constructing a global IMPES pressure equation from the MIIM equation, wherein said constructing 
the global IMPES pressure equation comprises: 

constructing a second scalar IMPES pressure equation for each of the implicit cells from die 
set of implicit equations corresponding to the implicit cell; and 
[0 concatenating the first scalar IMPES pressure equation fen* each of the IMPES cells and the 

second scalar IMPES pressure equation for each of the implicit cells; 

b) solving the global IMPES pressure equation for first pressures; 

c) computing improved saturations at the iix^)licit cells; 

d) determining if a convergence condition is satisfied; 

15 e) repeatedly perfonning b) through d) until the convergence condition is satisfied; 

f) computing a final solution estimate for die MIIM equation using the improved saturations and first 
pressures after the convergence condition is satisfied; 

g) applying the final solution estinutte to determine behavior of die reservoir model at a future discrete 
time value. 

20 

6. A method for performing reservoir simulation by solving an implicit linear equation arising in 
a Newton iteration of an implicit reservoir model, wherein the reservoir model comprises a plurality of cells, die 
method comprising: 

a) constructing a global IMPES pressure equation from the implicit Imear equation, wherein die 
25 global IMPES pressure equation comprises one scaUff IMPES pressure equation for each of die phuality of 

cells; 

b) solving the global IMPES pressure equation to determine first pressure values, wherein one of die 
first pressure values is associated with each of the pliirality of cells; 

c) constructing a complementary matrix equation in terms of unknowns other than pressure, wherein 
30 the complementary matrix equation is constructed using a constmint of conserving total velocity between cells; 

d) solving the complementary matrix equation to determine improved estimates of the unknowns 
odier than pressure at each of the plurality of cells; 

e) constructing a composite solution change which combines a first solution change associated with 
the first pressure values and a second solution change associated widi die improved estimates of unknowns other 

35 than pressure; 

f) providing die composite solution change to a solution accelerator, 

g) the solution accelerator generating an accelerated solution change; 

h) determining ifa convergence condition is satisfied; 

i) repeating (b) dirough (h) until the convergence condition is satisfied; 

40 j) coirqiuting a final solution estimate based on the accelerated solution change after the convergence 

condition is satisfiec^ 
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k) applying the final solution estinoate to predict properties of reservoir fluids at a future time value. 
7. The method of claim 6, wherein the solution accelerator is GMRES. 
5 8. The method ofclaott 6, wherein the sohition accelerator is ORTHOMIN. 

9. The method of claim 6, wherein said constructing a complementary matrix equation in terms 
of unknown other dian pressure comprises constructing a saturation matrix equation, wherein the unknowns 
other than pressure are saturations. 

10 

10. The method of claim 6, wherein said unkno\wns other dian pressure conqirise one or more 
saturations, mole fractions^ energies, masses, or volumes. 

11. The method of claim 6, further comprising computing first residuals of the implicit matrix 
15 equation based on the first pressure values, wherein the first residuals are used to construct the complementary 

matrix equatioa 

12. The method of claim 1 1, further comprising computing second residuals of the implicit matrix 
equation based on the improved estimates of die unknowns other than pressure, wherein the first residuals and 

20 second residuals are provided to die solution accelerator as input data. 

13. The method of claim 12, further comprising computing third residuals of die implicit matrix 
equation based on the accelerated solution change, wherein said determining if the convergence condition is 
satisfied con^rises determining if a magnitude of the diird residuals interpreted as a vector is snutUer than a 

25 threshold value. 

14. A method for performing reservoir simulation using total velocity sequential preconditioning, 
wherem the reservoir is sub-divided mto a phirality of cells, the method comprising: 

formulating finite difference equations which describe a behavior of reservoir fluids over a timestep; 
30 solving die finite difference equations by performing one or more Newton iterations, where each of 

said one or more Newton iteration conqprises: 

a) constructing a linear approximation for each non-linear term in the finite difference equations; 

b) constructing an iiiQ>licit matrix equation based on the finite difteience equations and die linear 
approxhnations; 

35 c) solving the implicit matrix equation, wherein said solving die implicit matrix equation com- 

prises: 

(ci) constructing a complementary matrix equation in terms of unknowns other than pressure 
usmg a constraint of conserving total velocity between cells; 

(c2) solving die complementary matrix equation for inqnoved estimates of unknowns other 

40 than pressure; 
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rqjeatedly pcrforaiing said solving the finite difference equations in order to predict behavior of the 
reservoir fluids over time. 

15. A mediod for perfonning reservoir simulation by solving an inq)lidt matrix equation arising 
5 from a Newton iteration of an inqslicit reservoir model, wherein die reservoir model con^rises a phirality of 

cells, wherein die in^licit matrix equation includes unknown variables, the method comprising: 

a) constructing a global IMPES pressure equation using die implicit matrix equation, wherein die 
global IMPES pressure equation comprises one scalar IMPES pressure equation for each of the plurality of 
cells; 

10 b) solving the global IMPES pressure equation to dcteimine first pressure values, wherein one of the 

first pressure vahies is associated with each of the plurality of cells; 

c) cQwi piting improved estimates of die unknown variables odier dian pressure by performing one or 
more iterations of a preconditionei; 

d) constructing a composite solution change by combining a first solution change associated with die 
15 first pressure values and a second solution change associated with die in^roved estimates of the unknown vari- 
ables odier than pressure; 

e) providing the composite solution change to a solution accelerator; 

Q the solution accelerator generating an accelerated solution change in response to the composite 
solution change; 

20 g) repeatedly perfonning b) through f) until a convergence criteria is satisfied; 

h) computing a final solution estimate based on die accelerated solution change after die conve^ence 
criteria is satisfied; 

wherein the final solution change is utilized to predict a behavior of reservoir fluids at a future discrete 
ixms value. 

25 

16. The method of claim 15, wherein die solution accelerator con^rises OrthomirL 

17. The mediod of claim 15, wherem die solution accelerator con^rises GMRES. 

30 18. A method for performing reservoir simulation by solving an implicit matrix equation arising 

from a Newton iteration of an inqjUcit reservoir model, wherein die implicit reservoir model comprises a plural- 
ity of ceUs, wherein die in^licit matrix equation in expressed m terms of unknown variables includmg pres- 
sures, the mediod conqirising: 

a) constructing a global IMPES pressure equation using die inq>licit matrix equation, wherein die 
35 global IMPES pressure equation comprises one scalar IMPES pressure equation for each of die phirality of 

cells; 

b) solving die global IMPES pressure equation to determine first improved estimates of die pressures, 
wherein one of die first improved estimates is associated widi each of the plurality of cells; 

c) confuting second inqmved estimates of all die unknowns variables by performing one or more 
40 iterations of a preconditioned 
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d) constructing a coniposite solution change in the unknoMoi variables by combining a first solution 
change associated with flie first improved estimates and a second solution change associated with the second 
in^roved estimatBs; 

e) providing die conqiosite sohition change to a solution accelerator; 

f) the solution accelerator generating an accelerated solution change 

g) determining if a convergence condition is satisfied; 

h) repeatedly performing b) through g) until the convergence condition is satisfied; 

i) computing a final solution estimate based on the accelerated solution change after the convergence 
condition is satisfied, and applying the final solution estimate to predict properties of reservoir fluids at a future 
time value. 

19. The method of claim 18, further comprising computing second residuals of the implicit matrix 
equation based on flie improved estimates of the unknown variables, wherein the second residuals are provided 
to the solution accelerator as input data. 

20. The method of claim 18, further comprising computing third residuals of die implicit matrix 
equation based on the accelerated solution change, wherein said determining if the convergence condition is 
satisfied con^rises determining if a magnitude of die third residual mterpreted as a vector is smaller dum a 
direshold value. 

21 . The method of claim 1 8, wherein the solution accelerator comprises OrthomixL 

22. The method of claim 18, wherein the solution accelerator comprises GMRES. 
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